---
title: "2023_regression_Electral voilence"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
library(haven)
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
#library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)

theme_set(theme_sjplot())

dat <- read_dta("data/Iraq2023_ver1.dta")


dat <- dat %>% 
  mutate(Distrust_gov = (q7_1_1 + q7_1_2 + q7_1_3 + q7_1_7) / 4,
         trust_gov = (trust_president + trust_PM + trust_parliament + trust_party) / 4) %>% 
  mutate(Distrust_gov2 = case_when(Distrust_gov <= 3 ~ "low",
                                   Distrust_gov >= 3.25 ~ "high"),
         trust_gov2 = case_when(trust_gov <= 3 ~ 0,
                                trust_gov >= 3.25 ~ 1)) %>%
  mutate(Minorities = others)

```

# overview
```{r}
stargazer(as.data.frame(dat), type = "text")
```

# overview dependent variavles
```{r}
hist(dat$violent_experiment)
dat %>% 
  filter(!is.na(violent_experiment)) %>% 
  ggplot() +
  geom_histogram(aes(x = violent_experiment, color = "white", bins = 4)) +
  labs(x = "Experiment", y = "frequency")
```

# regression: ethnicity
```{r}
reg01 <- lm(formula = violent_experiment ~ list_violent + sex, data = dat)
summary(reg01)

reg02 <- lm(formula = violent_experiment ~ list_violent*shia, data = dat)
summary(reg02)

reg03 <- lm(formula = violent_experiment ~ list_violent*sunni, data = dat)
summary(reg03)

reg04 <- lm(formula = violent_experiment ~ list_violent*kurd, data = dat)
summary(reg04)

reg05 <- lm(formula = violent_experiment ~ list_violent*others +
              sex + age + education + income, data = dat)
summary(reg05)

reg06 <- lm(formula = violent_experiment ~ list_violent*shia +
              sex + age + education + income, data = dat)
summary(reg06)

reg007 <- lm_robust(formula = violent_experiment ~ list_violent*shia +
                      sex + age + education + income, 
                   cluster = as.factor(C5), se_type = "stata",
                   data = dat)
summary(reg007) 
```

# regression ideology
```{r}
reg10 <- lm(formula = violent_experiment ~ list_violent*iraq_nationalism, data = dat)
summary(reg10)

reg11 <- lm(formula = violent_experiment ~ list_violent*iraq_nationalism +
              sex + age + education + income, data = dat)
summary(reg11)
```

# regression policy
```{r}
reg20 <- lm(formula = violent_experiment ~ list_violent*centralization, data = dat)
summary(reg20)

reg21 <- lm(formula = violent_experiment ~ list_violent*centralization +
              sex + age + education + income, data = dat)
summary(reg21)


reg22 <- lm(formula = violent_experiment ~ list_violent*independent_kurd, data = dat)
summary(reg22)

reg23 <- lm(formula = violent_experiment ~ list_violent*independent_kurd +
              sex + age + education + income, data = dat)
summary(reg23)

reg24 <- lm(formula = violent_experiment ~ list_violent * reconciliation +
              sex + age + education + income, data = dat)
summary(reg24)

reg25 <- lm(formula = violent_experiment ~ list_violent * anti_corruption +
              sex + age + education + income, data = dat)
summary(reg25)

```

# regression party
```{r}
reg30 <- lm(formula = violent_experiment ~ list_violent*Azm, data = dat)
summary(reg30)

reg31 <- lm(formula = violent_experiment ~ list_violent*Azm +
              sex + age + education + income, data = dat)
summary(reg31)

reg32 <- lm(formula = violent_experiment ~ list_violent*support_sunniparties +
               sex + age + education + income, data = dat)
summary(reg32)

reg33 <- lm(formula = violent_experiment ~ list_violent * support_governingshia,
            data = dat)
summary(reg33)

reg34 <- lm(formula = violent_experiment ~ list_violent * support_governingshia +
              sex + age + education + income,
            data = dat)
summary(reg34)

reg35 <- lm(formula = violent_experiment ~ list_violent * support_PMU_parties +
              sex + age + education + income,
            data = dat)
summary(reg35)

reg36 <- lm(formula = violent_experiment ~ list_violent * vote_governingshia +
              sex + age + education + income,
            data = dat)
summary(reg36)

reg37 <- lm(formula = violent_experiment ~ list_violent * vote_PMU_parties +
              sex + age + education + income,
            data = dat)
summary(reg37)
```

# regression trust
```{r}
reg40 <- lm(formula = violent_experiment ~ list_violent*trust_parliament, data = dat)
summary(reg40)

reg41 <- lm(formula = violent_experiment ~ list_violent*trust_parliament +
              sex + age + education + income, data = dat)
summary(reg41)

reg42 <- lm(formula = violent_experiment ~ list_violent*trust_religious_leader, data = dat)
summary(reg42)

reg42 <- lm(formula = violent_experiment ~ list_violent*trust_religious_leader +
            sex + age + education + income, data = dat)
summary(reg42)

reg43 <- lm(formula = violent_experiment ~ list_violent * distrust_gov +
            sex + age + education + income, data = dat)
summary(reg43)

reg44 <- lm(formula = violent_experiment ~ list_violent * Distrust_gov +
            sex + age + education + income, data = dat)
summary(reg44)

reg45 <- lm(formula = violent_experiment ~ list_violent * Distrust_gov2 +
            sex + age + education + income, data = dat)
summary(reg45)

reg46 <- lm(formula = violent_experiment ~ list_violent * trust_gov +
            sex + age + education + income, data = dat)
summary(reg46)

reg47 <- lm(formula = violent_experiment ~ list_violent * trust_gov2 +
            sex + age + education + income, data = dat)
summary(reg47)

```

# regression rely/ income
```{r}
reg50 <- lm(formula = violent_experiment ~ list_violent*rely_PMU_security, data = dat)
summary(reg50)

reg51 <- lm(formula = violent_experiment ~ list_violent*rely_PMU_robbery, data = dat)
summary(reg51)

reg52 <- lm(formula = violent_experiment ~ list_violent*rely_PMU_job, data = dat)
summary(reg52)

reg53 <- lm(formula = violent_experiment ~ list_violent*income, data = dat)
summary(reg53)

reg54 <- lm(formula = violent_experiment ~ list_violent*income +
              sex + age + education, data = dat)
summary(reg54)
```


# regression table
```{r}
reg_table <- mtable("M1(Treatment)" = reg01, "M2(Sect)" = reg05, "M3(Party)" = reg35, 
                    "M4(Trust)" = reg47, "M5(Policy)" = reg24, "M6(Policy)" = reg25,
               summary.stats=c("adj. R-squared","F","p", "N"))
print(reg_table, type = "html")
write_html(reg_table, file = "result/regression_table.html")
write.mtable(reg_table, file = "result/regression_table.csv", colsep = ",")
```

# coef_plot
```{r}
multiplot(reg01, reg05, reg35, reg47, reg24, reg25,
          title = "OLS results",
          names = c("M1:treatment", "M2:Sect", "M3:Party", 
                    "M4:Trust", "M5:Policy", "M6:Policy"),
          sort = c("natural"),
          intercept = FALSE, plot.linetypes = FALSE, legend.reverse = FALSE,
          plot.shapes = FALSE,
          numberAngle = 0, lwdOuter = 1, 
          coefficients = c("list_violent", 
                           "list_violent:others",
                           "list_violent:support_PMU_parties",
                           "list_violent:trust_gov2low",
                           "list_violent:reconciliation",
                           "list_violent:anti_corruption"),
          newNames = c(list_violent = "Treatment group", 
                       "list_violent:others" = "Minotiries",
                       "list_violent:support_PMU_parties" = "Support PMU parties",
                       "list_violent:trust_gov2low" = "Low trust of governments",
                       "list_violent:reconciliation" = "Reconciliation",
                       "list_violent:anti_corruption" = "Anti corruption"),
          single = FALSE,
          ncol = 3) +
  theme_bw() +
  theme(legend.position = "none")+
  ggtitle("")
```

# marginal effect: sect
```{r}
int_reg01 <- interplot(m = reg05, var1 = "list_violent", var2 = "others", point = TRUE) +
  labs(x = "Minorities",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(Model 1)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg01)

p1 <- plot_model(reg05, type = "pred", terms = c("others", "list_violent"), 
           axis.labels = "Treatment", title = "Subs. eff.(Model 1)",
           axis.title = c("Minorities", "Treatment"))
```

# marginal effect: party
```{r}
int_reg02 <- interplot(m = reg35, var1 = "list_violent", var2 = "support_PMU_parties", point = TRUE) +
  labs(x = "Support for PMU parties",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(Model 2)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg02)

p2 <- plot_model(reg35, type = "pred", terms = c("support_PMU_parties", "list_violent"), 
           axis.labels = "Treatment", title = "Subs. eff.(Model 2)",
           axis.title = c("Support for PMU parties", "Treatment"))
```

# marginal effect: trust
```{r}
int_reg03 <- interplot(m = reg47, var1 = "list_violent", var2 = "trust_gov2", point = TRUE) +
  labs(x = "Trust government",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(Model 3)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg03)

p3 <- plot_model(reg47, type = "pred", terms = c("trust_gov2", "list_violent"), 
           axis.labels = "Treatment", title = "Subs. eff.(Model 3)",
           axis.title = c("Trust government", "Treatment"))
```

# marginal effect: policy
```{r}
int_reg04 <- interplot(m = reg24, var1 = "list_violent", var2 = "reconciliation", point = TRUE) +
  labs(x = "Reconciliation",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(Model 4)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg04)

p4 <- plot_model(reg24, type = "pred", terms = c("reconciliation", "list_violent"), 
           axis.labels = "Treatment", title = "Subs. eff.(Model 4)",
           axis.title = c("Reconciliation", "Treatment"))
```

# marginal effect: policy
```{r}
int_reg05 <- interplot(m = reg25, var1 = "list_violent", var2 = "anti_corruption", point = TRUE) +
  labs(x = "Anti corruption",
       y = "Effect of treatment") +
  ggtitle("Marg. eff.(Model 5)") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg05)

p5 <- plot_model(reg25, type = "pred", terms = c("anti_corruption", "list_violent"), 
           axis.labels = "Treatment", title = "Subs. eff.(Model 5)",
           axis.title = c("Anti corruption", "Treatment"))
```

# plot substantial effect
```{r}
mf <- ggarrange(int_reg01, 
          int_reg02 + theme(axis.title.y = element_blank()), 
          int_reg03 + theme(axis.title.y = element_blank()),
          int_reg04, 
          int_reg05 + theme(axis.title.y = element_blank()), 
          nrow = 2, ncol = 3, align = "hv")
mf

ggsave(filename = "result/marginal_effect.png",
       plot = mf,
       width = 8,
       height = 5,
       units = "in", 
       dpi = 600,
       bg = "white")
```
# plot marginal effect
```{r}
sf <- ggarrange(p1, 
          p2 + theme(axis.title.y = element_blank()), 
          p3 + theme(axis.title.y = element_blank()),
          p4, 
          p5 +theme(axis.title.y = element_blank()),
          nrow = 2, ncol = 3, align = "hv",
          common.legend = TRUE,
          legend = "right")
sf

ggsave(filename = "result/substantial_effect.png",
       plot = sf,
       width = 8,
       height = 5,
       units = "in", 
       dpi = 600,
       bg = "white")
```


# robustness check
```{r}
reg200 <- lm(formula = violent_experiment ~ list_violent + sex, data = dat)
summary(reg200)

reg201 <- lm(formula = violent_experiment ~ list_violent, 
             data = subset(dat, sex == 0))
summary(reg201)

reg202 <- lm(formula = violent_experiment ~ list_violent, 
             data = subset(dat, sex == 1))
summary(reg202)

# Inverse Probability Weighting
ps_model <- glm(list_violent ~ sex,
                data = dat,
                family = binomial())
dat$pscore <- predict(ps_model, type = "response")

dat <- dat |>
  mutate(weight_ipw = ifelse(list_violent == 1, 1 / pscore, 1/(1 - pscore)),
         treat_mean = mean(list_violent),
         weight_stab = ifelse(list_violent == 1, treat_mean / pscore, (1 - treat_mean)/ (1 - pscore)))

model_ipw <- lm(violent_experiment ~ list_violent, 
                data = dat,
                weights = weight_ipw)
summary(model_ipw)

model_stab <- lm(violent_experiment ~ list_violent, 
                data = dat,
                weights = weight_stab)
summary(model_stab)

# table
robust_table <- mtable("Controled" = reg200, "Male subset" = reg201, "Female subset" = reg202, 
                       "IPW" = model_ipw, "Stability" = model_stab,
               summary.stats=c("adj. R-squared","F","p", "N"))
print(robust_table, type = "html")
write_html(robust_table, file = "result/robust_table.html")
write.mtable(robust_table, file = "result/robust_table.csv", colsep = ",")
```


# region
```{r}
reg100 <- lm(formula = violent_experiment ~ list_violent * North +
              sex + age + education, data = dat)
summary(reg100)

reg101 <- lm(formula = violent_experiment ~ list_violent * West +
              sex + age + education, data = dat)
summary(reg101)

reg102 <- lm(formula = violent_experiment ~ list_violent * Central +
              sex + age + education, data = dat)
summary(reg102)

reg103 <- lm(formula = violent_experiment ~ list_violent * South +
              sex + age + education, data = dat)
summary(reg103)

reg104 <- lm(formula = violent_experiment ~ list_violent * KRG +
              sex + age + education, data = dat)
summary(reg104)


reg_region <- mtable("North" = reg100, "West" = reg101, 
                     "Central" = reg102, "South" = reg103, "KRG" = reg104,
                     summary.stats = c("adj. R-squared","F","p", "N"))
print(reg_region, type = "html")
write_html(reg_region, file = "data/result_reigon2023.html")
```

# marginal effect: region
```{r}
int_reg07 <- interplot(m = reg100, var1 = "list_violent", var2 = "North", point = TRUE) +
  labs(x = "North***",
       y = "Effect of treatment") +
  ggtitle("Model 6") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg07)

int_reg08 <- interplot(m = reg101, var1 = "list_violent", var2 = "West", point = TRUE) +
  labs(x = "West***",
       y = "Effect of treatment") +
  ggtitle("Model 7") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg08)

int_reg09 <- interplot(m = reg102, var1 = "list_violent", var2 = "Central", point = TRUE) +
  labs(x = "Central",
       y = "Effect of treatment") +
  ggtitle("Model 8") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg09)

int_reg10 <- interplot(m = reg103, var1 = "list_violent", var2 = "South", point = TRUE) +
  labs(x = "South***",
       y = "Effect of treatment") +
  ggtitle("Model 9") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg10)

int_reg11 <- interplot(m = reg104, var1 = "list_violent", var2 = "KRG", point = TRUE) +
  labs(x = "KRG***",
       y = "Effect of treatment") +
  ggtitle("Model 10") +
  geom_hline(yintercept = 0, linetype = "dashed")
print(int_reg11)

ggarrange(int_reg07, 
          int_reg08 + theme(axis.title.y = element_blank()), 
          int_reg09 + theme(axis.title.y = element_blank()), 
          int_reg10, 
          int_reg11 + theme(axis.title.y = element_blank()),
          nrow = 2, ncol = 3, align = "hv")
```
